Covariance2D Subroutine

private subroutine Covariance2D(dist, theta, cov, sill, range, varmodel)

Subroutine to calculate a covariance matrix of observations from a semivariogram model and a vector of separation distances.

Arguments

Type IntentOptional Attributes Name
real(kind=float), intent(in), DIMENSION(:,:) :: dist

separation distance vector
real(kind=float), intent(in), DIMENSION(:,:) :: theta

anisotropy angle ??
real(kind=float), intent(out), DIMENSION(:,:) :: cov

covariance matrix
real(kind=float), intent(in) :: sill

partial sill (total sill - nugget) from automatic fitting
real(kind=float), intent(in) :: range

range from automatic fitting
integer(kind=short), intent(in) :: varmodel

semivariogram model: 1 = spherical, 2 = exponential, 3 = gaussian, 4 = power

Variables

Type Visibility Attributes Name Initial
real(kind=float), public, DIMENSION(SIZE(dist,1),SIZE(dist,2)) :: h_prime
integer(kind=short), public :: i

Source Code

SUBROUTINE Covariance2D &
!
(dist, theta, cov, sill, range, varmodel)

IMPLICIT NONE

!Arguments with intent(in)
REAL (KIND = float), INTENT(in) ::  sill  !! partial sill (total sill - nugget) from automatic fitting
REAL (KIND = float), INTENT(in) :: range  !! range from automatic fitting
REAL (KIND = float), DIMENSION(:,:), INTENT(IN) :: dist !!separation distance vector
REAL (KIND = float), DIMENSION(:,:), INTENT(IN) :: theta !!anisotropy angle ??
INTEGER (KIND = short), INTENT(in) :: varmodel  !!semivariogram model: 1 = spherical, 2 = exponential, 3 = gaussian, 4 = power

!Arguments with intent(out):
REAL (KIND = float), DIMENSION(:,:), INTENT(OUT) :: cov  !!covariance matrix

!Local variable declarations
REAL (KIND = float), DIMENSION(SIZE(dist,1),SIZE(dist,2)) :: h_prime 
INTEGER (KIND = short) :: i

!---------------------------------end of declarations--------------------------

!h_prime = dist ! da correggere con anisotropia quando sarà implementata

	
!Check for anisotropy
IF (anisotropyAngle == 0.) THEN
	h_prime = dist
ELSE
	!Transform distance vector by applying anisotropy correction factor
 	h_prime = dist + dist * ABS( SIN( theta - anisotropyAngle ) * anisotropyRatio )
END IF

!covariance = nugget + partialsill - semivariogram(h)
!since semivariogram(h) = nugget + partialsill * f(h,range) => covariance = partialsill * ( 1- f(h,range) )
SELECT CASE(varmodel)
	CASE (1) !Spherical
		WHERE(h_prime < range)
        cov = sill * (1. - ( 1.5 * (h_prime / range) - 0.5 * (h_prime / range)**3. ) )
		ELSEWHERE
        cov = 0.
        END WHERE
             
    CASE (2) !Exponential
        WHERE(h_prime < range)
        !cov = sill * (1. - ( 1. - EXP (- h_prime / range) ) )
        cov = sill * (1. - ( 1. - EXP (- 3*  h_prime / range) ) )
		ELSEWHERE
        cov = 0.
        END WHERE
             
    CASE (3) !Gaussian: MUST have a nugget effect to function without instability
        WHERE(h_prime < range)
        !cov = sill * (1. - ( 1. - EXP (- h_prime * h_prime / (range * range) ) ) ) 
        cov = sill * (1. - ( 1. - EXP (- 10 * h_prime * h_prime / (range * range) ) ) ) 
		ELSEWHERE
        cov = 0.
        END WHERE
             
    CASE (4) !Power
        WHERE(h_prime < range)
        cov = sill * (1. - (h_prime / range) ) 
		ELSEWHERE
        cov = 0.
        END WHERE

END SELECT 

RETURN	
END SUBROUTINE Covariance2D